Análisis exploratorio univariante

Se tiene dataset consitutido por observaciones de 11 variables en 34 estados del mundo. Las variables se encuentran estandarizadas pues están tomadas con unidades de medida muy diferentes. Las variables son:

Mostramos los datos:

library(foreign)
data <- read.spss("./DB_3.sav", to.data.frame = TRUE)
data

Comprobamos que las variables están normalizadas:

summary(data)
##      PAIS              ZPOBDENS          ZTMINFAN          ZESPVIDA      
##  Length:34          Min.   :-1.0778   Min.   :-1.1026   Min.   :-2.1453  
##  Class :character   1st Qu.:-0.8497   1st Qu.:-0.9586   1st Qu.:-0.7460  
##  Mode  :character   Median :-0.1616   Median :-0.3931   Median : 0.2781  
##                     Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000  
##                     3rd Qu.: 0.5547   3rd Qu.: 0.8985   3rd Qu.: 0.9033  
##                     Max.   : 2.8616   Max.   : 1.9048   Max.   : 1.2486  
##                                                                          
##     ZPOBURB           ZTMEDICO          ZPAGRICU          ZPSERVI        
##  Min.   :-1.7697   Min.   :-1.1473   Min.   :-1.2342   Min.   :-1.88521  
##  1st Qu.:-0.7320   1st Qu.:-0.8829   1st Qu.:-0.8480   1st Qu.:-0.72858  
##  Median : 0.1268   Median :-0.2916   Median :-0.2134   Median : 0.03541  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.00000  
##  3rd Qu.: 0.8148   3rd Qu.: 0.8694   3rd Qu.: 0.7115   3rd Qu.: 0.85176  
##  Max.   : 1.5096   Max.   : 2.3717   Max.   : 1.9052   Max.   : 1.62885  
##                                                                          
##     ZTLIBROP          ZTEJERCI           ZTPOBACT          ZTENERGI      
##  Min.   :-0.9696   Min.   :-0.86586   Min.   :-2.1341   Min.   :-0.9507  
##  1st Qu.:-0.9240   1st Qu.:-0.59889   1st Qu.:-0.6735   1st Qu.:-0.7813  
##  Median :-0.3237   Median :-0.20626   Median :-0.1067   Median :-0.3900  
##  Mean   : 0.0000   Mean   : 0.00000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.: 0.7736   3rd Qu.: 0.07996   3rd Qu.: 0.8789   3rd Qu.: 0.5828  
##  Max.   : 2.4024   Max.   : 4.42620   Max.   : 1.7045   Max.   : 2.7498  
##  NA's   :1

Mostramos desviación típica de las variables

print("Desviación típica:")
## [1] "Desviación típica:"
sapply(data, sd)
## Warning in var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm =
## na.rm): NAs introduced by coercion
##     PAIS ZPOBDENS ZTMINFAN ZESPVIDA  ZPOBURB ZTMEDICO ZPAGRICU  ZPSERVI 
##       NA        1        1        1        1        1        1        1 
## ZTLIBROP ZTEJERCI ZTPOBACT ZTENERGI 
##       NA        1        1        1

Desviación típica de la variable ZTLIBROP sin el valor perdido:

print("Desviación típica:")
## [1] "Desviación típica:"
sd(data[,9][-22])
## [1] 1

Recodificaciones o agrupaciones

Recodificamos a UTF-8 y aprovechamos para corregir el nombre de algunos paises incompletos y/o con erratas.

data[1, "PAIS"] <- "áfrica sur"
data[3, "PAIS"] <- "argentina"
data[4, "PAIS"] <- "australia"
data[11, "PAIS"]<- "españa"
data[14, "PAIS"]<- "hungría"
data[16, "PAIS"]<- "indonesia"
data[20, "PAIS"]<- "japón"
data[22, "PAIS"]<- "marruecos"
data[23, "PAIS"] <- "méxico"
data[27, "PAIS"] <- "rd alemania"
data[28, "PAIS"] <- "reino unido"
data[29, "PAIS"] <- "rf alemania"

Valores perdidos

Vemos los valores perdidos:

cbind(apply(is.na(data),2,sum),apply(is.na(data),2,sum)/dim(data)[1])
##          [,1]       [,2]
## PAIS        0 0.00000000
## ZPOBDENS    0 0.00000000
## ZTMINFAN    0 0.00000000
## ZESPVIDA    0 0.00000000
## ZPOBURB     0 0.00000000
## ZTMEDICO    0 0.00000000
## ZPAGRICU    0 0.00000000
## ZPSERVI     0 0.00000000
## ZTLIBROP    1 0.02941176
## ZTEJERCI    0 0.00000000
## ZTPOBACT    0 0.00000000
## ZTENERGI    0 0.00000000

vemos un valor perdido en ZTLIBROP. Aprovechamos para guardar una copia del dataset sin variable PAIS que no es de interés.

data_ <- data[-1]
head(data_)

Los datos perdidos en ZTLIBROP son inferiores al 5%, sustituiremos por la media:

not_available<-function(data,na.rm=F){
  data[is.na(data)]<-mean(data,na.rm=T)
  data
}

data_<-as.data.frame(apply(data_, 2, not_available))

Análisis descriptivo numérico clásico

Veamos un resumen de las variables:

summary(data_)
##     ZPOBDENS          ZTMINFAN          ZESPVIDA          ZPOBURB       
##  Min.   :-1.0778   Min.   :-1.1026   Min.   :-2.1453   Min.   :-1.7697  
##  1st Qu.:-0.8497   1st Qu.:-0.9586   1st Qu.:-0.7460   1st Qu.:-0.7320  
##  Median :-0.1616   Median :-0.3931   Median : 0.2781   Median : 0.1268  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.: 0.5547   3rd Qu.: 0.8985   3rd Qu.: 0.9033   3rd Qu.: 0.8148  
##  Max.   : 2.8616   Max.   : 1.9048   Max.   : 1.2486   Max.   : 1.5096  
##     ZTMEDICO          ZPAGRICU          ZPSERVI            ZTLIBROP      
##  Min.   :-1.1473   Min.   :-1.2342   Min.   :-1.88521   Min.   :-0.9696  
##  1st Qu.:-0.8829   1st Qu.:-0.8480   1st Qu.:-0.72858   1st Qu.:-0.9145  
##  Median :-0.2916   Median :-0.2134   Median : 0.03541   Median :-0.2442  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.00000   Mean   : 0.0000  
##  3rd Qu.: 0.8694   3rd Qu.: 0.7115   3rd Qu.: 0.85176   3rd Qu.: 0.6923  
##  Max.   : 2.3717   Max.   : 1.9052   Max.   : 1.62885   Max.   : 2.4024  
##     ZTEJERCI           ZTPOBACT          ZTENERGI      
##  Min.   :-0.86586   Min.   :-2.1341   Min.   :-0.9507  
##  1st Qu.:-0.59889   1st Qu.:-0.6735   1st Qu.:-0.7813  
##  Median :-0.20626   Median :-0.1067   Median :-0.3900  
##  Mean   : 0.00000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.: 0.07996   3rd Qu.: 0.8789   3rd Qu.: 0.5828  
##  Max.   : 4.42620   Max.   : 1.7045   Max.   : 2.7498

Mostramos boxplot de las variables para visualizar su distribución:

boxplot (data_, main = "Boxplot de las variables", xlab = "Variables", ylab = "Valores")

Mostramos coef. de asimetría de las variables

library(moments)
skewness(data_)
##   ZPOBDENS   ZTMINFAN   ZESPVIDA    ZPOBURB   ZTMEDICO   ZPAGRICU    ZPSERVI 
##  1.0679539  0.5477565 -0.5802673 -0.3189850  0.4909934  0.5918253 -0.1239223 
##   ZTLIBROP   ZTEJERCI   ZTPOBACT   ZTENERGI 
##  0.8170779  2.8653904 -0.1284517  1.2270338

Mostramos curtosis de las variables

kurtosis(data_)
##  ZPOBDENS  ZTMINFAN  ZESPVIDA   ZPOBURB  ZTMEDICO  ZPAGRICU   ZPSERVI  ZTLIBROP 
##  3.542159  1.796001  2.156273  1.902481  2.065480  2.055664  1.911288  2.524152 
##  ZTEJERCI  ZTPOBACT  ZTENERGI 
## 12.634918  2.111551  3.912472

Analizaremos con más detalle:

ZPOBDENS

## $centralidad
## $centralidad$clasica
## $centralidad$clasica$media
## [1] 1.06373e-16
## 
## 
## $centralidad$resistente
##                   [,1]
## PMC         -0.1474835
## Trimedia    -0.1545405
## Centrimedia -0.2293829
## 
## 
## $dispersion
## $dispersion$clasica
## $dispersion$clasica$desviación_típica
## [1] 1
## 
## $dispersion$clasica$Coef_varización
## [1] 9.400883e+15
## 
## $dispersion$clasica$rango
## [1] -1.077834  2.861621
## 
## 
## $dispersion$resistente
##                               75%
## Rango Inter-Cuartílico  1.4043955
## MEDA                    0.6924864
## CVc                    -4.7611940
## 
## 
## $forma
## $forma$clasica
## $forma$clasica$skewness
## [1] 1.067954
## 
## $forma$clasica$kurtosis
## [1] 3.542159
## 
## 
## $forma$resistente
##                                         25%
## Asimetría de Yule               -0.08733974
## Asimetría de Kelly              -0.37369403
## Asimetría de Kelly adimensional  2.31250000

Las medidas resistentes de centralidad están ligeramente hacia la izquierda. MEDA inferior a desviación típica. Observamos asimetría, así como cierta acumulación de datos por el valor de curtosis.

ZTMINFAN

## $centralidad
## $centralidad$clasica
## $centralidad$clasica$media
## [1] 8.982296e-17
## 
## 
## $centralidad$resistente
##                   [,1]
## PMC         -0.0300511
## Trimedia    -0.2115862
## Centrimedia -0.2268480
## 
## 
## $dispersion
## $dispersion$clasica
## $dispersion$clasica$desviación_típica
## [1] 1
## 
## $dispersion$clasica$Coef_varización
## [1] 1.113301e+16
## 
## $dispersion$clasica$rango
## [1] -1.102554  1.904824
## 
## 
## $dispersion$resistente
##                                75%
## Rango Inter-Cuartílico   1.8571199
## MEDA                     0.6040459
## CVc                    -30.8993711
## 
## 
## $forma
## $forma$clasica
## $forma$clasica$skewness
## [1] 0.5477565
## 
## $forma$clasica$kurtosis
## [1] 1.796001
## 
## 
## $forma$resistente
##                                        25%
## Asimetría de Yule               -0.9235577
## Asimetría de Kelly              -0.6132994
## Asimetría de Kelly adimensional  1.5600769

Observamos cierto desplazamiento a la izquierda en las medidas resistentes de centralidad. Asimetría algo menor que el caso anterior, así como una distribución más aplanada.

ZESPVIDA

## $centralidad
## $centralidad$clasica
## $centralidad$clasica$media
## [1] -1.72884e-16
## 
## 
## $centralidad$resistente
##                   [,1]
## PMC         0.07863105
## Trimedia    0.17836465
## Centrimedia 0.17613181
## 
## 
## $dispersion
## $dispersion$clasica
## $dispersion$clasica$desviación_típica
## [1] 1
## 
## $dispersion$clasica$Coef_varización
## [1] -5.784225e+15
## 
## $dispersion$clasica$rango
## [1] -2.145279  1.248640
## 
## 
## $dispersion$resistente
##                               75%
## Rango Inter-Cuartílico  1.6493257
## MEDA                    0.7621433
## CVc                    10.4877506
## 
## 
## $forma
## $forma$clasica
## $forma$clasica$skewness
## [1] -0.5802673
## 
## $forma$clasica$kurtosis
## [1] 2.156273
## 
## 
## $forma$resistente
##                                        25%
## Asimetría de Yule               -0.7172544
## Asimetría de Kelly               0.4995611
## Asimetría de Kelly adimensional  1.7963476

Observamos distrbución aplanada aunque menos que la anterior, variación en los datos, así como medidas resistentes de centralidad ligeramente desplazadas a la derecha.

ZPOBURB

## $centralidad
## $centralidad$clasica
## $centralidad$clasica$media
## [1] 4.511206e-16
## 
## 
## $centralidad$resistente
##                  [,1]
## PMC         0.0413792
## Trimedia    0.0840905
## Centrimedia 0.1147624
## 
## 
## $dispersion
## $dispersion$clasica
## $dispersion$clasica$desviación_típica
## [1] 1
## 
## $dispersion$clasica$Coef_varización
## [1] 2.216702e+15
## 
## $dispersion$clasica$rango
## [1] -1.769694  1.509616
## 
## 
## $dispersion$resistente
##                               75%
## Rango Inter-Cuartílico  1.5467796
## MEDA                    0.7223655
## CVc                    18.6903015
## 
## 
## $forma
## $forma$clasica
## $forma$clasica$skewness
## [1] -0.318985
## 
## $forma$clasica$kurtosis
## [1] 1.902481
## 
## 
## $forma$resistente
##                                        25%
## Asimetría de Yule               -0.6736702
## Asimetría de Kelly               0.2958259
## Asimetría de Kelly adimensional  2.3329787

Observamos distribución con cierta asimetría, la más aplanada tras la de ZTMINFAN, y prácticamente no hay desplazamiento de centralidad.

ZTMEDICO

## $centralidad
## $centralidad$clasica
## $centralidad$clasica$media
## [1] 2.416495e-17
## 
## 
## $centralidad$resistente
##                     [,1]
## PMC         -0.006716126
## Trimedia    -0.149133349
## Centrimedia -0.149734266
## 
## 
## $dispersion
## $dispersion$clasica
## $dispersion$clasica$desviación_típica
## [1] 1
## 
## $dispersion$clasica$Coef_varización
## [1] 4.138224e+16
## 
## $dispersion$clasica$rango
## [1] -1.147256  2.371712
## 
## 
## $dispersion$resistente
##                                 75%
## Rango Inter-Cuartílico    1.7522727
## MEDA                      0.7980172
## CVc                    -130.4526316
## 
## 
## $forma
## $forma$clasica
## $forma$clasica$skewness
## [1] 0.4909934
## 
## $forma$clasica$kurtosis
## [1] 2.06548
## 
## 
## $forma$resistente
##                                        25%
## Asimetría de Yule               -0.9769641
## Asimetría de Kelly              -0.3735297
## Asimetría de Kelly adimensional  1.2811833

Se aprecia distribución aplanada, cierta asimetría, así como ligero desplazamiento a la izquierda.

ZPAGRICU

## $centralidad
## $centralidad$clasica
## $centralidad$clasica$media
## [1] 6.355338e-17
## 
## 
## $centralidad$resistente
##                    [,1]
## PMC         -0.06825485
## Trimedia    -0.14081091
## Centrimedia -0.20433276
## 
## 
## $dispersion
## $dispersion$clasica
## $dispersion$clasica$desviación_típica
## [1] 1
## 
## $dispersion$clasica$Coef_varización
## [1] 1.57348e+16
## 
## $dispersion$clasica$rango
## [1] -1.234234  1.905157
## 
## 
## $dispersion$resistente
##                                75%
## Rango Inter-Cuartílico   1.5595319
## MEDA                     0.7069276
## CVc                    -11.4243309
## 
## 
## $forma
## $forma$clasica
## $forma$clasica$skewness
## [1] 0.5918253
## 
## $forma$clasica$kurtosis
## [1] 2.055664
## 
## 
## $forma$resistente
##                                        25%
## Asimetría de Yule               -0.6801059
## Asimetría de Kelly              -0.4817497
## Asimetría de Kelly adimensional  2.2578456

Distribución aplanada, cierta asimetría, así como medidas resistentes de centralidad ligeramente desplazadas a la izquierda.

ZPSERVI

## $centralidad
## $centralidad$clasica
## $centralidad$clasica$media
## [1] -3.529078e-17
## 
## 
## $centralidad$resistente
##                     [,1]
## PMC         0.0615893232
## Trimedia    0.0485015920
## Centrimedia 0.0006495749
## 
## 
## $dispersion
## $dispersion$clasica
## $dispersion$clasica$desviación_típica
## [1] 1
## 
## $dispersion$clasica$Coef_varización
## [1] -2.833601e+16
## 
## $dispersion$clasica$rango
## [1] -1.885211  1.628845
## 
## 
## $dispersion$resistente
##                               75%
## Rango Inter-Cuartílico  1.5803435
## MEDA                    0.7918077
## CVc                    12.8296875
## 
## 
## $forma
## $forma$clasica
## $forma$clasica$skewness
## [1] -0.1239223
## 
## $forma$clasica$kurtosis
## [1] 1.911288
## 
## 
## $forma$resistente
##                                        25%
## Asimetría de Yule               0.73913043
## Asimetría de Kelly              0.05071496
## Asimetría de Kelly adimensional 1.43206522

Distribución aplanada, prácticamente centrada, y poca asimetría.

ZTLIBROP

## $centralidad
## $centralidad$clasica
## $centralidad$clasica$media
## [1] 6.729069e-17
## 
## 
## $centralidad$resistente
##                   [,1]
## PMC         -0.1111009
## Trimedia    -0.1776580
## Centrimedia -0.2615358
## 
## 
## $dispersion
## $dispersion$clasica
## $dispersion$clasica$desviación_típica
## [1] 0.9847319
## 
## $dispersion$clasica$Coef_varización
## [1] 1.4634e+16
## 
## $dispersion$clasica$rango
## [1] -0.9696156  2.4023604
## 
## 
## $dispersion$resistente
##                               75%
## Rango Inter-Cuartílico  1.6067690
## MEDA                    0.6849277
## CVc                    -7.2311230
## 
## 
## $forma
## $forma$clasica
## $forma$clasica$skewness
## [1] 0.8170779
## 
## $forma$clasica$kurtosis
## [1] 2.524152
## 
## 
## $forma$resistente
##                                        25%
## Asimetría de Yule               -0.5450695
## Asimetría de Kelly              -0.5339910
## Asimetría de Kelly adimensional  2.1865595

Se observa fuerte asimetría hacia la derecha, así como desplazamiento a la derecha y concentración de los datos por el alto valor de curtosis. ZTEJERCI

## $centralidad
## $centralidad$clasica
## $centralidad$clasica$media
## [1] 1.859184e-16
## 
## 
## $centralidad$resistente
##                   [,1]
## PMC         -0.2594621
## Trimedia    -0.2328606
## Centrimedia -0.2192510
## 
## 
## $dispersion
## $dispersion$clasica
## $dispersion$clasica$desviación_típica
## [1] 1
## 
## $dispersion$clasica$Coef_varización
## [1] 5.378705e+15
## 
## $dispersion$clasica$rango
## [1] -0.8658606  4.4262018
## 
## 
## $dispersion$resistente
##                               75%
## Rango Inter-Cuartílico  0.6788462
## MEDA                    0.3313997
## CVc                    -1.3081800
## 
## 
## $forma
## $forma$clasica
## $forma$clasica$skewness
## [1] 2.86539
## 
## $forma$clasica$kurtosis
## [1] 12.63492
## 
## 
## $forma$resistente
##                                        25%
## Asimetría de Yule                0.2579423
## Asimetría de Kelly              -0.3448882
## Asimetría de Kelly adimensional  1.6721112

Observamos muy fuerte asímetría a la derecha, desplazamiento de las medidas resistentes de centralidad a la izquierda, así como alta concentración de los datos por valor muy alto de curtosis. ZTPOBACT

## $centralidad
## $centralidad$clasica
## $centralidad$clasica$media
## [1] -2.952027e-16
## 
## 
## $centralidad$resistente
##                    [,1]
## PMC          0.10268877
## Trimedia    -0.00198850
## Centrimedia  0.02917055
## 
## 
## $dispersion
## $dispersion$clasica
## $dispersion$clasica$desviación_típica
## [1] 1
## 
## $dispersion$clasica$Coef_varización
## [1] -3.387502e+15
## 
## $dispersion$clasica$rango
## [1] -2.134094  1.704472
## 
## 
## $dispersion$resistente
##                              75%
## Rango Inter-Cuartílico 1.5523378
## MEDA                   0.8557514
## CVc                    7.5584589
## 
## 
## $forma
## $forma$clasica
## $forma$clasica$skewness
## [1] -0.1284517
## 
## $forma$clasica$kurtosis
## [1] 2.111551
## 
## 
## $forma$resistente
##                                         25%
## Asimetría de Yule               -1.96271531
## Asimetría de Kelly              -0.06440751
## Asimetría de Kelly adimensional  0.60382547

Distribución aplanada, poca asimetría y prácticamente centrada. ZTENERGI

## $centralidad
## $centralidad$clasica
## $centralidad$clasica$media
## [1] 3.586955e-17
## 
## 
## $centralidad$resistente
##                    [,1]
## PMC         -0.09924855
## Trimedia    -0.24461072
## Centrimedia -0.26271176
## 
## 
## $dispersion
## $dispersion$clasica
## $dispersion$clasica$desviación_típica
## [1] 1
## 
## $dispersion$clasica$Coef_varización
## [1] 2.78788e+16
## 
## $dispersion$clasica$rango
## [1] -0.950661  2.749785
## 
## 
## $dispersion$resistente
##                              75%
## Rango Inter-Cuartílico  1.364062
## MEDA                    0.520457
## CVc                    -6.871948
## 
## 
## $forma
## $forma$clasica
## $forma$clasica$skewness
## [1] 1.227034
## 
## $forma$clasica$kurtosis
## [1] 3.912472
## 
## 
## $forma$resistente
##                                        25%
## Asimetría de Yule               -0.7454988
## Asimetría de Kelly              -0.5056537
## Asimetría de Kelly adimensional  1.2966382

Se observa distribución algo desplazada a la izquierda, concentrada (alto valor de curtosis), y asimétrica a la derecha.

Valores extremos (outliers)

Veamos de nuevo boxplot de las variables:

boxplot (data_, main = "Boxplot de las variables", xlab = "Variables", ylab = "Valores")

En cierta manera todas las variables parecen tener cierta similitud a una normal. Vemos 1 outlier en ZPOBDENS y ZTENERGI, varios en ZTEJERCI. Los sustituiremos por la media.

outlier<-function(data,na.rm=T){
  H<-1.5*IQR(data)
  
  if(any(data<=(quantile(data,0.25,na.rm = T)-H))){
    data[data<=quantile(data,0.25,na.rm = T)-H]<-NA
    data[is.na(data)]<-mean(data,na.rm=T)
    data<-outlier(data)}
  
  if(any(data>=(quantile(data,0.75, na.rm = T)+H))){
    data[data>=quantile(data,0.75, na.rm = T)+H]<-NA
    data[is.na(data)]<-mean(data,na.rm=T)
    data<-outlier(data)
  }
  return(data)
}

data_<-apply(data_, 2, outlier)
boxplot(data_, main = "Boxplot de las variables", xlab = "Variables", ylab = "Valores")

Mostramos resumen de los datos tras tratamiento de outliers:

summary(data_)
##     ZPOBDENS           ZTMINFAN          ZESPVIDA          ZPOBURB       
##  Min.   :-1.07783   Min.   :-1.1026   Min.   :-2.1453   Min.   :-1.7697  
##  1st Qu.:-0.84968   1st Qu.:-0.9586   1st Qu.:-0.7460   1st Qu.:-0.7320  
##  Median :-0.16160   Median :-0.3931   Median : 0.2781   Median : 0.1268  
##  Mean   :-0.08672   Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.: 0.40244   3rd Qu.: 0.8985   3rd Qu.: 0.9033   3rd Qu.: 0.8148  
##  Max.   : 2.15204   Max.   : 1.9048   Max.   : 1.2486   Max.   : 1.5096  
##     ZTMEDICO          ZPAGRICU          ZPSERVI            ZTLIBROP      
##  Min.   :-1.1473   Min.   :-1.2342   Min.   :-1.88521   Min.   :-0.9696  
##  1st Qu.:-0.8829   1st Qu.:-0.8480   1st Qu.:-0.72858   1st Qu.:-0.9145  
##  Median :-0.2916   Median :-0.2134   Median : 0.03541   Median :-0.2442  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.00000   Mean   : 0.0000  
##  3rd Qu.: 0.8694   3rd Qu.: 0.7115   3rd Qu.: 0.85176   3rd Qu.: 0.6923  
##  Max.   : 2.3717   Max.   : 1.9052   Max.   : 1.62885   Max.   : 2.4024  
##     ZTEJERCI           ZTPOBACT          ZTENERGI      
##  Min.   :-0.86586   Min.   :-2.1341   Min.   :-0.9507  
##  1st Qu.:-0.59889   1st Qu.:-0.6735   1st Qu.:-0.7813  
##  Median :-0.28969   Median :-0.1067   Median :-0.3900  
##  Mean   :-0.28969   Mean   : 0.0000   Mean   :-0.1690  
##  3rd Qu.:-0.04621   3rd Qu.: 0.8789   3rd Qu.: 0.3346  
##  Max.   : 0.68708   Max.   : 1.7045   Max.   : 1.5672

y los datos iniciales:

summary(data)
##      PAIS              ZPOBDENS          ZTMINFAN          ZESPVIDA      
##  Length:34          Min.   :-1.0778   Min.   :-1.1026   Min.   :-2.1453  
##  Class :character   1st Qu.:-0.8497   1st Qu.:-0.9586   1st Qu.:-0.7460  
##  Mode  :character   Median :-0.1616   Median :-0.3931   Median : 0.2781  
##                     Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000  
##                     3rd Qu.: 0.5547   3rd Qu.: 0.8985   3rd Qu.: 0.9033  
##                     Max.   : 2.8616   Max.   : 1.9048   Max.   : 1.2486  
##                                                                          
##     ZPOBURB           ZTMEDICO          ZPAGRICU          ZPSERVI        
##  Min.   :-1.7697   Min.   :-1.1473   Min.   :-1.2342   Min.   :-1.88521  
##  1st Qu.:-0.7320   1st Qu.:-0.8829   1st Qu.:-0.8480   1st Qu.:-0.72858  
##  Median : 0.1268   Median :-0.2916   Median :-0.2134   Median : 0.03541  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.00000  
##  3rd Qu.: 0.8148   3rd Qu.: 0.8694   3rd Qu.: 0.7115   3rd Qu.: 0.85176  
##  Max.   : 1.5096   Max.   : 2.3717   Max.   : 1.9052   Max.   : 1.62885  
##                                                                          
##     ZTLIBROP          ZTEJERCI           ZTPOBACT          ZTENERGI      
##  Min.   :-0.9696   Min.   :-0.86586   Min.   :-2.1341   Min.   :-0.9507  
##  1st Qu.:-0.9240   1st Qu.:-0.59889   1st Qu.:-0.6735   1st Qu.:-0.7813  
##  Median :-0.3237   Median :-0.20626   Median :-0.1067   Median :-0.3900  
##  Mean   : 0.0000   Mean   : 0.00000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.: 0.7736   3rd Qu.: 0.07996   3rd Qu.: 0.8789   3rd Qu.: 0.5828  
##  Max.   : 2.4024   Max.   : 4.42620   Max.   : 1.7045   Max.   : 2.7498  
##  NA's   :1

Supuesto de normalidad

Veamos si estamos tratando con variables normales, para ello utilizamos el método gráfico qq-plot.

par(mar=c(1,1,1,1))
par(mfrow=c(3,5))
invisible(apply(data_, 2, function(x){
  qqnorm(x,main=NULL)
  abline(a=0,b=1,col="red")
}))

Vemos que las variables que más se acercan a la normalidad son la 4, 6, 7 y 10 (ZPOBURB, ZPAGRICU, ZPSERVI y ZTPOBACT), en menor medida la 1, 2, 3, 5, 8 y 11 (ZPOBDENS, ZTMINFAN, ZESPVIDA, ZTMEDICO, ZTLIBROP y ZTENERGI), la que menos se acerca a normalidad es la 9 (ZTEJERCI).

Normalizamos los datos

data_norm <- scale(data_)

Análisis exploratorio multivariante

Trabajaremos con los datos normalizados:

datos_pca <- data_norm
row.names(datos_pca) = data$PAIS

Correlación entre variables, test de Barlett

Mostramos las correlaciones entre variables

library(corrr)
datos_pca_correlaciones <- correlate(datos_pca)  #Cálculo de objeto de correlaciones
## 
## Correlation method: 'pearson'
## Missing treated using: 'pairwise.complete.obs'
rplot(datos_pca_correlaciones, legend = TRUE, colours = c("firebrick1", "black","darkcyan"), print_cor = TRUE) 
## Don't know how to automatically pick scale for object of type noquote. Defaulting to continuous.
## Don't know how to automatically pick scale for object of type noquote. Defaulting to continuous.

## 
## Attaching package: 'psych'
## The following object is masked _by_ '.GlobalEnv':
## 
##     outlier

Hacemos test de esfericidad de Bartlett por el que comprobaremos si las correlaciones son significativamente distintas de 0

cortest.bartlett(cor(datos_pca), nrow(datos_pca))
## $chisq
## [1] 387.5628
## 
## $p.value
## [1] 1.583951e-51
## 
## $df
## [1] 55

Observamos un p-valor muy bajo, se rechaza hipótesis nula por lo que los datos no están incorrelados.

Haremos Análisis de Componentes Principales, los datos ya los hemos preparado para ello.

PCA<-prcomp(datos_pca)

Mostramos los coeficientes de las componentes principales, esto es, el peso de cada variable en la correspondiente componente principal:

PCA$rotation
##                  PC1         PC2         PC3          PC4         PC5
## ZPOBDENS  0.06796131  0.39036079  0.07578571 -0.873056137  0.05981546
## ZTMINFAN -0.37529798 -0.11573027 -0.08212636 -0.014741875  0.41788865
## ZESPVIDA  0.37165513  0.07270402 -0.00407411  0.012970861 -0.53281544
## ZPOBURB   0.35808785 -0.29124050 -0.06094609 -0.113795410 -0.05152603
## ZTMEDICO  0.33374417  0.13909259 -0.09587119  0.280930962  0.15782496
## ZPAGRICU -0.36224935  0.29940891  0.04348437  0.061235615 -0.15998056
## ZPSERVI   0.29607632 -0.48878963  0.09183567 -0.148420697 -0.03064997
## ZTLIBROP  0.32594252  0.08892433 -0.07078859 -0.135187927  0.28062939
## ZTEJERCI  0.02397325  0.17948476 -0.93319975 -0.006241743 -0.12518016
## ZTPOBACT  0.20289865  0.56167354  0.29659750  0.281932179 -0.11796713
## ZTENERGI  0.33155119  0.20152127 -0.02054058  0.148110844  0.61275026
##                  PC6         PC7         PC8          PC9        PC10
## ZPOBDENS -0.24337351  0.09619809 -0.04353317  0.022846323 -0.01646372
## ZTMINFAN -0.01756878  0.05245475 -0.39097670 -0.157260993 -0.47091956
## ZESPVIDA  0.03028763 -0.04105430  0.29289664 -0.091799724 -0.54354628
## ZPOBURB  -0.07065324 -0.10651593 -0.29868358 -0.792450091  0.16228759
## ZTMEDICO -0.36060911  0.78403566 -0.06560314  0.001929436 -0.06125296
## ZPAGRICU  0.13248858  0.09343760 -0.01200167 -0.325592955 -0.48343682
## ZPSERVI  -0.07352269 -0.09280074 -0.37422057  0.449245272 -0.39418278
## ZTLIBROP  0.85707761  0.21913126 -0.01229527  0.015266047 -0.02755826
## ZTEJERCI -0.02870986 -0.16141551 -0.18218984  0.127931321  0.01085377
## ZTPOBACT  0.02948798 -0.25103160 -0.61023870  0.097328148  0.08409493
## ZTENERGI -0.21278600 -0.45367225  0.34643517 -0.066593066 -0.23535405
##                  PC11
## ZPOBDENS -0.016686681
## ZTMINFAN -0.511610948
## ZESPVIDA -0.424556791
## ZPOBURB   0.088144761
## ZTMEDICO  0.058614024
## ZPAGRICU  0.617855535
## ZPSERVI   0.361609762
## ZTLIBROP  0.008370356
## ZTEJERCI  0.063899183
## ZTPOBACT -0.095710436
## ZTENERGI  0.143880098

Mostramos las desviaciones típicas de cada componente principal, proporción de varianza explicada y acumulada:

PCA$sdev
##  [1] 2.4818030 1.2805818 1.0311152 0.9515420 0.6448055 0.6027181 0.4818699
##  [8] 0.3354192 0.2491626 0.1641974 0.1390789
summary(PCA)
## Importance of components:
##                           PC1    PC2     PC3     PC4    PC5     PC6     PC7
## Standard deviation     2.4818 1.2806 1.03112 0.95154 0.6448 0.60272 0.48187
## Proportion of Variance 0.5599 0.1491 0.09665 0.08231 0.0378 0.03302 0.02111
## Cumulative Proportion  0.5599 0.7090 0.80568 0.88799 0.9258 0.95881 0.97992
##                            PC8     PC9    PC10    PC11
## Standard deviation     0.33542 0.24916 0.16420 0.13908
## Proportion of Variance 0.01023 0.00564 0.00245 0.00176
## Cumulative Proportion  0.99015 0.99579 0.99824 1.00000

Mostramos gráficamente la varianza explicada

## 
## Attaching package: 'ggplot2'
## The following objects are masked from 'package:psych':
## 
##     %+%, alpha
varianza_explicada <- PCA$sdev^2 / sum(PCA$sdev^2)
ggplot(data = data.frame(varianza_explicada, pc = 1:11),
       aes(x = pc, y = varianza_explicada, fill=varianza_explicada )) +
  geom_col(width = 0.3) +
  scale_y_continuous(limits = c(0,0.6)) + theme_bw() +
  labs(x = "Componente principal", y= " Proporción de varianza explicada")

Y la varianza acumulada

varianza_acum<-cumsum(varianza_explicada)
ggplot( data = data.frame(varianza_acum, pc = 1:11),
        aes(x = pc, y = varianza_acum ,fill=varianza_acum )) +
  geom_col(width = 0.5) +
  scale_y_continuous(limits = c(0,1)) +
  theme_bw() +
  labs(x = "Componentes principales",
       y = "Proporción varianza acumulada")

Reducción de dimensión mediante variables observables

Para la selección de componentes principales utilizaremos la regla de Abdi et al. (2010). Promediamos las varianzas explicadas por la componentes principales:

PCA$sdev^2
##  [1] 6.15934602 1.63988972 1.06319852 0.90543213 0.41577420 0.36326910
##  [7] 0.23219856 0.11250601 0.06208202 0.02696079 0.01934294

Calculamos la media:

mean(PCA$sdev^2)
## [1] 1

Y no quedamos con las tres primeras componentes principales pues tienen varianza explicada superior a 1, que es la media.

Visualizaremos ahora la contribución de las componentes principales.

Veamos una comparativa entre la primera y segunda componente principal analizando que variables tienen más peso para la definición de cada componente principal

library("factoextra")
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
fviz_pca_var(PCA,
             repel=TRUE,col.var="cos2",
             legend.title="Distancia")+theme_bw()

Ahora comparativa entre la primera y tercera componente principal analizando que variables tienen más peso para la definición de cada componente principal

fviz_pca_var(PCA,axes=c(1,3),
             repel=TRUE,col.var="cos2",
             legend.title="Distancia")+theme_bw()

Comparativa entre la segunda y tercera componente principal analizando que variables tienen más peso para la definición de cada componente principal

fviz_pca_var(PCA,axes=c(2,3),
             repel=TRUE,col.var="cos2",
             legend.title="Distancia")+theme_bw()

Representaremos las observaciones de los objetos junto con las componentes principales, identificando con colores las observaciones que mayor varianza explican de las componentes principales.

Observaciones en la primera y segunda componente principal

fviz_pca_ind(PCA,col.ind = "contrib",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel=TRUE,legend.title="Contrib.var")+theme_bw()

Observaciones en la primera y tercera componente principal

fviz_pca_ind(PCA,axes=c(1,3),col.ind = "contrib",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel=TRUE,legend.title="Contrib.var")+theme_bw()

Observaciones en la segunda y tercera componente principal

fviz_pca_ind(PCA,axes=c(2,3),col.ind = "contrib",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel=TRUE,legend.title="Contrib.var")+theme_bw()

Veremos ahora representaciones conjuntas de individuos y variables, en las que se muestran las contribuciones de los individuos a las componentes principales y el peso de cada variable en las componentes principales.

Variables y observaciones en la primera y segunda componente principal

fviz_pca(PCA,
         alpha.ind ="contrib", col.var = "cos2",col.ind="seagreen",
         gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"),
         repel=TRUE,
         legend.title="Distancia")+theme_bw()

Variables y observaciones en la primera y tercera componente principal

fviz_pca(PCA,axes=c(1,3),
         alpha.ind ="contrib", col.var = "cos2",col.ind="seagreen",
         gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"),
         repel=TRUE,
         legend.title="Distancia")+theme_bw()

Variables y observaciones en la segunda y tercera componente principal

fviz_pca(PCA,axes=c(2,3),
         alpha.ind ="contrib", col.var = "cos2",col.ind="seagreen",
         gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"),
         repel=TRUE,
         legend.title="Distancia")+theme_bw()

Mostramos las coordenadas de los datos originales tipificados en el nuevo sistema de referencia.

PCA$x
##                    PC1         PC2         PC3         PC4         PC5
## áfrica sur  -1.1879510 -1.57174870  0.45044485  0.37814927  1.42283283
## argelia     -2.3462462 -1.88106620 -0.81081919  0.20145035  0.45280276
## argentina    1.1979236 -1.53833675  0.20077839  0.89727380 -0.38850796
## australia    3.2068509 -0.89415852  0.83940147  1.06527240  0.55028348
## brasil      -0.8909634 -1.34735175  0.89893299  0.55508246 -0.13891942
## canada       2.7566836 -1.13584507  1.63380445  0.83558450 -0.40111280
## chile        0.3748013 -2.00974285 -0.30168062  0.13027061 -1.15642746
## china       -2.1745978  2.61356042  0.96542176  0.79447745 -1.51186129
## coreasur    -0.2395061  0.22577834 -0.02200379 -0.29949146 -0.46326227
## egipto      -2.8332447 -0.89472527 -1.74869384 -0.02886629  0.06798311
## españa       1.9771834 -0.23696615 -1.33663589  0.04438475 -0.47649830
## filipina    -2.1452831  0.20351988  1.27403374 -0.86565659 -0.56838811
## francia      2.6414111 -0.01330041 -0.62271201  0.03170101 -0.08010665
## hungría      1.8369158  1.71813286 -0.93757556  0.47001001  0.46001652
## india       -3.7665089  1.25869208  1.25543964 -1.00446829  0.51111537
## indonesia   -3.5122382 -0.02386264  1.26366487  0.10062250  0.42626703
## iran        -2.1106915 -0.95775804 -0.92365774  0.22851570 -0.20529497
## israel       2.5185220 -1.13346621 -0.22444498 -0.43374508 -0.19339430
## italia       1.8793472  0.33682331 -0.08587639 -0.53435523 -0.44626390
## japón        2.4042017  0.77519469  1.67063103 -2.16353425 -0.37124266
## libano       0.3908113 -0.93862865 -1.80755715 -2.54094256 -0.39063809
## marruecos   -2.6111083 -0.63494215 -0.84370677 -0.04560012  0.32314717
## méxico      -0.4292694 -1.88248069  1.06058617  0.18898306 -0.18870059
## nigeria     -4.0691890  0.51312721  1.11817945  0.17003627  0.74676184
## pakistan    -3.7752500 -0.29164848 -0.45582169 -0.43747966  0.81514836
## polonia      1.4051850  1.79382646 -0.34557393  0.76060587  0.06958040
## rd alemania  3.0006719  1.80429843 -0.54847315  0.41154507  0.79574301
## reino unido  3.3366283  0.19074891  0.77266251 -1.61255080  0.44377012
## rf alemania  3.6113690  0.93736385 -0.47800299 -1.52413606  1.06149382
## rumania      0.7068599  1.72531630 -0.57683789  0.83431076  0.06419030
## turquia     -2.2589795  1.18954704 -2.39282243  0.64594436 -0.65460455
## urss         2.0656447  0.85194470  0.06708247  2.12196985  0.75589667
## usa          2.1968825 -0.80816808  0.55718509  0.86676472 -0.76556968
## vietnam     -3.1568661  2.05632216  0.43464714 -0.24212838 -0.56623978
##                     PC6          PC7         PC8         PC9        PC10
## áfrica sur  -0.27199158 -0.623595597 -0.09144981  0.34753493  0.36296233
## argelia     -0.21767407 -0.113018266  0.67613362  0.31890906 -0.22213535
## argentina   -0.77144470  0.807228821 -0.15129705 -0.40998101  0.13438946
## australia    0.38518213 -0.811803953  0.25022260 -0.16193690 -0.38067934
## brasil       0.02175312  0.175309634 -0.23014686 -0.55123723  0.07607427
## canada       1.26624990  0.214600937 -0.41474660  0.39540041 -0.09652726
## chile       -0.18275397 -0.754674607 -0.13156616 -0.14262580  0.15878588
## china        0.21360710 -0.315617302  0.16287419  0.09008596 -0.04968074
## coreasur     1.11284552 -0.035491221  0.39806886 -0.23092496  0.26947888
## egipto       0.02250445 -0.310866567 -0.08354177 -0.02901969 -0.19432866
## españa       0.82678614  0.777742199  0.48167714 -0.10099465  0.01181882
## filipina    -0.15028934  0.004493093  0.49597736  0.08712967  0.12559592
## francia      0.13854922 -0.260637363  0.10989111  0.31541253 -0.16761311
## hungría      0.88238111  0.833538559  0.01574888  0.25974597  0.17058796
## india       -0.20212478  0.579046678 -0.03383508 -0.25659043 -0.21257778
## indonesia    0.13318440  0.115415733 -0.06815793  0.33036118 -0.08581324
## iran         0.12404861 -0.318863831  0.50032561 -0.17821369  0.21184100
## israel       0.10225569  0.802635989 -0.30878270 -0.14148837 -0.06317868
## italia      -1.29071596  0.829905483  0.37519670  0.17010251 -0.10240096
## japón       -0.68185453 -0.375862213 -0.00357082  0.01804949 -0.07676834
## libano      -0.69319566  0.064675345 -0.38259552  0.27888287  0.13462981
## marruecos    0.86093321  0.077375097  0.12651960 -0.09999321 -0.03204079
## méxico      -0.56822183 -0.038134015  0.30275495 -0.16956075 -0.16756015
## nigeria      0.22834537  0.182754443 -0.44626412 -0.22970715  0.15560703
## pakistan    -0.13809459  0.282765202 -0.24412230  0.23118350 -0.06426054
## polonia     -0.49281049 -0.275363943  0.12995559 -0.09238346  0.01584053
## rd alemania -0.84486979 -0.720954096 -0.03279271 -0.22607825  0.14754470
## reino unido  0.60953258 -0.592216125 -0.31921611 -0.18569892 -0.00228887
## rf alemania  0.72445007  0.104964263  0.11641759 -0.19173155 -0.05768685
## rumania     -0.20141636 -0.406726190  0.36289963  0.11225930  0.12367151
## turquia      0.03651040 -0.430409682 -0.82424809 -0.19586549 -0.26359945
## urss        -0.97941596  0.614086633 -0.12622830  0.03093806 -0.01039851
## usa         -0.07839738 -0.061084782 -0.51133189  0.47419412  0.12224111
## vietnam      0.07615199 -0.021218358 -0.10076961  0.13384194  0.02846942
##                      PC11
## áfrica sur   0.0009993838
## argelia     -0.1452388044
## argentina   -0.0562911808
## australia    0.2406669888
## brasil      -0.2189263852
## canada      -0.1070909000
## chile        0.1209991593
## china        0.0175549674
## coreasur     0.1594077183
## egipto      -0.0573909279
## españa       0.1049016670
## filipina     0.2772966071
## francia      0.0703157933
## hungría     -0.0509025203
## india       -0.1058064933
## indonesia    0.0175770992
## iran        -0.2075043694
## israel       0.0268765029
## italia      -0.0920630334
## japón       -0.2402706122
## libano       0.2033940439
## marruecos   -0.2011568416
## méxico       0.1298877117
## nigeria      0.1771866008
## pakistan     0.0721410202
## polonia      0.0187492876
## rd alemania -0.0453764779
## reino unido -0.0604243295
## rf alemania  0.0292461854
## rumania     -0.1884743097
## turquia      0.0158359250
## urss         0.1909211317
## usa         -0.1407301850
## vietnam      0.0436895774

Reducción de dimensión mediante variables latentes

Como comprobamos anteriormente las correlaciones son significativamente distintas de 0, por lo que tiene sentido el Análisis Factorial.

Podemos visualizar las correlaciones:

## 
## Attaching package: 'polycor'
## The following object is masked from 'package:psych':
## 
##     polyserial
data_fa <- data_norm
poly_cor <- hetcor(data_fa)$correlations
ggcorrplot(poly_cor, type = "lower", hc.order = T)

Ahora, compararemos las salidas de modelos de 3 factores con el método del factor principal y con el de máxima verosimilitud.

modelo1<-fa(poly_cor,
            nfactors = 3,
            rotate = "none",
            fm="mle") # modelo máxima verosimilitud
modelo2<-fa(poly_cor,
            nfactors = 3,
            rotate = "none",
            fm="pa") # modelo factor principal
## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
## ultra-Heywood case was detected. Examine the results carefully

Ante la advertencia decidimos probar el método de mínimo residuo

modelo2<-fa(poly_cor,
            nfactors = 3,
            rotate = "none",
            fm="minres") # modelo mínimo residuo

Finalmente probamos otro método de mínimo residuo

modelo2<-fa(poly_cor,
            nfactors = 3,
            rotate = "none",
            fm="old.min") # modelo mínimo residuo
## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.

Comparamos las comunalidades:

sort(modelo1$communality,decreasing = T)->c1
sort(modelo2$communality,decreasing = T)->c2
head(cbind(c1,c2))
##                 c1        c2
## ZPAGRICU 0.9896670 0.9862799
## ZTMINFAN 0.9768999 0.9560850
## ZESPVIDA 0.9632148 0.9192056
## ZTENERGI 0.9577942 0.8752384
## ZPSERVI  0.9286098 0.8726959
## ZPOBURB  0.9120259 0.7414501

Comparamos las unicidades, es decir, la proporción de varianza que no ha sido explicada por el factor (1-comunalidad):

sort(modelo1$uniquenesses,decreasing = T)->u1
sort(modelo2$uniquenesses,decreasing = T)->u2
head(cbind(u1,u2))
##                  u1        u2
## ZTEJERCI 0.97397969 0.9449885
## ZPOBDENS 0.90544077 0.4806882
## ZTLIBROP 0.41684561 0.4225841
## ZTMEDICO 0.31612567 0.4110643
## ZTPOBACT 0.23641177 0.2962771
## ZPOBURB  0.08797406 0.2585499

Determinaremos el número óptimo de factores por método Scree plot (Cattel 1966) y Análisis Paralelo (Horn 1965):

scree(poly_cor)

fa.parallel(poly_cor, n.obs = nrow(data_fa), fa = "fa", fm = "mle")

## Parallel analysis suggests that the number of factors =  2  and the number of components =  NA

El número óptimo de factores es 2.

Estimamos el modelo factorial con 2 factores implementando una rotación tipo varimax para buscar una interpretación más simple.

modelo_varimax<-fa(poly_cor,nfactors = 2,rotate = "varimax", fa="mle")

Mostramos la matriz de pesos factorial rotada

print(modelo_varimax$loadings,cut=0) 
## 
## Loadings:
##          MR1    MR2   
## ZPOBDENS  0.023  0.285
## ZTMINFAN -0.756 -0.574
## ZESPVIDA  0.769  0.509
## ZPOBURB   0.964  0.057
## ZTMEDICO  0.631  0.524
## ZPAGRICU -0.985 -0.051
## ZPSERVI   0.941 -0.245
## ZTLIBROP  0.647  0.425
## ZTEJERCI  0.000  0.115
## ZTPOBACT  0.135  0.839
## ZTENERGI  0.603  0.575
## 
##                  MR1   MR2
## SS loadings    5.143 2.239
## Proportion Var 0.468 0.204
## Cumulative Var 0.468 0.671

Veamos una representación más interpretable

fa.diagram(modelo_varimax)

Vemos que el primer factor está asociado con ZPAGRICU, ZPOBURB,ZPSERVI, ZESPVIDA, ZTMINFAN, ZTLIBROP, ZTMEDICO y ZTENERGI, mientras que el segundo factor está asociado con ZTPOBACT.

Veamos con el test de hipótesis que contrasta si el número de factores es suficiente.

factanal(data_fa,factors=2, rotation="none")
## 
## Call:
## factanal(x = data_fa, factors = 2, rotation = "none")
## 
## Uniquenesses:
## ZPOBDENS ZTMINFAN ZESPVIDA  ZPOBURB ZTMEDICO ZPAGRICU  ZPSERVI ZTLIBROP 
##    0.908    0.020    0.053    0.071    0.402    0.051    0.085    0.446 
## ZTEJERCI ZTPOBACT ZTENERGI 
##    0.987    0.341    0.453 
## 
## Loadings:
##          Factor1 Factor2
## ZPOBDENS  0.147   0.266 
## ZTMINFAN -0.960  -0.244 
## ZESPVIDA  0.953   0.195 
## ZPOBURB   0.889  -0.373 
## ZTMEDICO  0.763   0.128 
## ZPAGRICU -0.887   0.403 
## ZPSERVI   0.758  -0.584 
## ZTLIBROP  0.743         
## ZTEJERCI          0.111 
## ZTPOBACT  0.460   0.669 
## ZTENERGI  0.732   0.104 
## 
##                Factor1 Factor2
## SS loadings      5.884   1.298
## Proportion Var   0.535   0.118
## Cumulative Var   0.535   0.653
## 
## Test of the hypothesis that 2 factors are sufficient.
## The chi square statistic is 63.09 on 34 degrees of freedom.
## The p-value is 0.00176

Vemos p valor muy bajo, se rechaza la hipótesis nula, 2 no son suficientes.

Estimamos el modelo factorial con 3 factores

Veamos una representación más interpretable

modelo_varimax<-fa(poly_cor,nfactors = 3,rotate = "varimax", fa="mle")
fa.diagram(modelo_varimax)

Vemos que el primer factor está asociado con ZPAGRICU, ZPOBURB,ZPSERVI, ZESPVIDA, ZTMINFAN, ZTLIBROP y ZTMEDICO, mientras que el segundo factor está asociado con ZTPOBACT y ZTENERGI, el tercero con ZTEJERCI.

Veamos con el test de hipótesis que contrasta si el número de factores es suficiente.

factanal(data_fa,factors=3, rotation="none")
## 
## Call:
## factanal(x = data_fa, factors = 3, rotation = "none")
## 
## Uniquenesses:
## ZPOBDENS ZTMINFAN ZESPVIDA  ZPOBURB ZTMEDICO ZPAGRICU  ZPSERVI ZTLIBROP 
##    0.905    0.023    0.037    0.088    0.316    0.010    0.071    0.417 
## ZTEJERCI ZTPOBACT ZTENERGI 
##    0.974    0.236    0.042 
## 
## Loadings:
##          Factor1 Factor2 Factor3
## ZPOBDENS          0.291         
## ZTMINFAN -0.893  -0.395   0.154 
## ZESPVIDA  0.888   0.354  -0.221 
## ZPOBURB   0.930  -0.199         
## ZTMEDICO  0.760   0.269   0.185 
## ZPAGRICU -0.961   0.255         
## ZPSERVI   0.831  -0.455  -0.178 
## ZTLIBROP  0.738   0.180         
## ZTEJERCI          0.119   0.105 
## ZTPOBACT  0.357   0.775   0.186 
## ZTENERGI  0.770   0.281   0.535 
## 
##                Factor1 Factor2 Factor3
## SS loadings      5.915   1.477   0.487
## Proportion Var   0.538   0.134   0.044
## Cumulative Var   0.538   0.672   0.716
## 
## Test of the hypothesis that 3 factors are sufficient.
## The chi square statistic is 26.21 on 25 degrees of freedom.
## The p-value is 0.397

No se rechaza la hipótesis nula de que sean suficientes.

Análisis de la normalidad multivariante

Utilizaremos funciones del paquete MVN para contrastar la normalidad multivariante. Veamos si tenemos outliers:

library(MVN)
outliers <- mvn(data = datos_pca  , mvnTest = "hz", multivariateOutlierMethod = "quan")

Se detectan 11 outliers en las observaciones. Veamos tests de Royston y Henze-Zirkler:

royston_test <- mvn(data = datos_pca , mvnTest = "royston", multivariatePlot = "qq")

royston_test$multivariateNormality
hz_test <- mvn(data = datos_pca , mvnTest = "hz")
hz_test$multivariateNormality

El test de Royston nos indica rechazar la hip. nula, no tendríamos normalidad multivariante. El test de Henze-Zirkler no nos sugiere rechzar la hip. nula de normalidad multivariante, aunque da p-valor bajo.

Clustering

# install.packages("tidyverse")
# install.packages("cluster")
# install.packages("factoextra")
# Cargamos los paquetes indicados
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ tibble  3.1.4     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.3     ✓ stringr 1.4.0
## ✓ readr   2.0.1     ✓ forcats 0.5.1
## ✓ purrr   0.3.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x ggplot2::%+%()   masks psych::%+%()
## x ggplot2::alpha() masks psych::alpha()
## x dplyr::filter()  masks stats::filter()
## x dplyr::lag()     masks stats::lag()
library(cluster)
library(factoextra)
data_clust <- datos_pca

Matriz de distancias:

distance<- get_dist(data_clust)
fviz_dist(distance, gradient = list(low ="#00AFBB", mid = "white", high = "#FC4E07"))

Aplicamos clustering con K-Medias:

k2 <- kmeans(data_clust, centers = 2, nstart = 25)
k3 <- kmeans(data_clust, centers = 3, nstart = 25)
k4 <- kmeans(data_clust, centers = 4, nstart = 25)
k5 <- kmeans(data_clust, centers = 5, nstart = 25)
k6 <- kmeans(data_clust, centers = 6, nstart = 25)

# plots to compare
p1 <- fviz_cluster(k2, geom = "point", data = data_clust) + ggtitle("k = 2")
p2 <- fviz_cluster(k3, geom = "point", data = data_clust) + ggtitle("k = 3")
p3 <- fviz_cluster(k4, geom = "point", data = data_clust) + ggtitle("k = 4")
p4 <- fviz_cluster(k5, geom = "point", data = data_clust) + ggtitle("k = 5")
p5 <- fviz_cluster(k6, geom = "point", data = data_clust) + ggtitle("k = 6")
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
grid.arrange(p1, p2, p3, p4, p5, nrow = 2)

Buscando el número óptimo de clusters:

set.seed(123)
fviz_nbclust(data_clust, kmeans, method = "wss")

Método de Silhouette

fviz_nbclust(data_clust, kmeans, method = "silhouette")

Método estadístico de brecha (GAP)

gap_stat <- clusGap(data_clust, FUN = kmeans, nstart = 25, K.max = 10, B = 50)
fviz_gap_stat(gap_stat)

Dos de los métodos sugieren que K=2 es el número de clusters más adecuado.

set.seed(123)
final <- kmeans(data_clust, 2, nstart = 25)
print(final)
## K-means clustering with 2 clusters of sizes 18, 16
## 
## Cluster means:
##     ZPOBDENS   ZTMINFAN   ZESPVIDA    ZPOBURB   ZTMEDICO   ZPAGRICU    ZPSERVI
## 1  0.2043031 -0.7853941  0.7769316  0.7240592  0.7671309 -0.7448802  0.5956415
## 2 -0.2298410  0.8835684 -0.8740481 -0.8145666 -0.8630222  0.8379902 -0.6700966
##     ZTLIBROP   ZTEJERCI   ZTPOBACT   ZTENERGI
## 1  0.6366902  0.1499519  0.4474332  0.6725565
## 2 -0.7162765 -0.1686959 -0.5033623 -0.7566261
## 
## Clustering vector:
##  áfrica sur    argelia    argentina   australia    brasil      canada   
##           2           2           1           1           2           1 
##    chile       china       coreasur    egipto        españa    filipina 
##           1           2           2           2           1           2 
##    francia      hungría    india      indonesia    iran        israel   
##           1           1           2           2           2           1 
##    italia         japón    libano     marruecos      méxico    nigeria  
##           1           1           1           2           2           2 
##    pakistan    polonia  rd alemania reino unido rf alemania    rumania  
##           2           1           1           1           1           1 
##    turquia     urss        usa         vietnam  
##           2           1           1           2 
## 
## Within cluster sum of squares by cluster:
## [1] 107.32705  88.71239
##  (between_SS / total_SS =  46.0 %)
## 
## Available components:
## 
## [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
## [6] "betweenss"    "size"         "iter"         "ifault"
# Visualizamos los resultados
fviz_cluster(final, data = data_clust)

Finalmente, mostramos medias de las variables a nivel de cluster.

as.data.frame(data_clust) %>%
  mutate(Cluster = final$cluster) %>%
  group_by(Cluster) %>%
  summarise_all("mean")